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Multi-particle dynamics in one-dimensional asymmetric exclusion processes with disorder is inves- 
tigated theoretically by computational and analytical methods. It is argued that the general phase 
diagram consists of three non-equilibrium phases that are determined by the dynamic behavior at 
the entrance, at the exit and at the slowest defect bond in the bulk of the system. Specifically, we 
consider dynamics of asymmetric exclusion process with two identical defect bonds as a function 
of distance between them. Two approximate theoretical methods, that treat the system as a se- 
quence of segments with exact description of dynamics inside the segments and neglect correlations 
between them, are presented. In addition, a numerical iterative procedure for calculating dynamic 
properties of asymmetric exclusion systems is developed. Our theoretical predictions are compared 
with extensive Monte Carlo computer simulations. It is shown that correlations play an important 
role in the particle dynamics. When two defect bonds are far away from each other the strongest 
correlations are found at these bonds. However, bringing defect bonds closer leads to the shift of 
correlations to the region between them. Our analysis indicates that it is possible to develop a suc- 
cessful theoretical description of asymmetric exclusion processes with disorder by properly taking 
into account the correlations. 
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I. INTRODUCTION 



In recent years a significant attention has been devoted 
to investigation of low-dimensional asymmetric simple 
exclusion processes (ASEPs) [J 0, SL S . They play a 
critical role for understanding fundamental properties of 
non-equilibrium phenomena in Chemistry, Physics and 
Biology. ASEPs have been widely utilized for descrip- 
tion of traffic phenomena [4L kinetics of biopolymcriza- 
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tion [6||, protein synthesis 
transport of motor proteins 
using asymmetric exclusion processes for studying mech- 
anisms of non-equilibrium phenomena is due to the fact 
that some homogeneous versions of ASEPs can be solved 
exactly via matrix-product approach and related meth- 
ods [ElaO- I n addition, understanding of processes in 
ASEPs can be achieved by utilizing a phenomenological 
domain wall approach [l4| . In order to have a more real- 
istic description of different non-equilibrium phenomena 
ASEPs with inhomogeneous distribution of rates are re- 
quired. However, there is a limited number of studies 
dealing with ASEPs with disorder in the transition rates 



at sites (static impurities) 



Mi m, m 



So, 11 HI M HE m, M MMMMmS^^ with 

disorder associated to particles 's hopping rates (moving 



* Corresponding author: foolad@iasbs.ac.ir 



impurities) [3J, [3f| [3y, [37|, [38[ • In this case exact solu- 
tions are not obtained, and extensive Monte Carlo com- 
puter simulations and approximate theories are utilized 
in order to understand particle dynamics. Disorder has 
a strong effect on the behavior of ASEPs. Even a single 
defect bond far away from the boundaries lead to dra- 
matic effects in the stationary properties both in closed 
[IH [l|| and open boundary conditions . It was shown 
recently that the dynamics of ASEPs is also influenced by 
several defects that are close to each other [IMS^al- 
though the mechanism of this phenomenon is not well un- 
derstood. This interaction between defects is important 
for understanding several biological transport phenom- 
ena [1, G3 ■ Recently the particular case of two defects 
has been extensively investigated by Monte Carlo simu- 
lations jlQ ]. It has been shown that the system current 
exhibits a notable dependence on the distance between 
defects with equal hopping rates. Moreover, it was found 
that the density profile is linear between defects which 
marks the existence of wandering shock between defects 
[n|. The case of two defective sites with equal rate has 
been generalized to include extended objects [3l|. Theo- 
retical efforts to analyze ASEPs with disorder have been 
m ostly directed to the cases with a single or few defects 
B Gil SI- In Ref. ASEP with open boundaries and 
with a local inhomogeneity in the bulk has been investi- 
gated by arguing that the defect bond divides the system 
into two coupled homogeneous ASEPs. This theoretical 
approach can be called a defect mean-field (DMF) be- 
cause the mean-field assumptions are made only at the 
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position of local inhomogeneity. Although a good agree- 
ment with computer simulations has been found, there 
were significant deviations in statistical properties of the 
phase with the maximal current that was attributed to 
the neglect of correlations at the defect bond in the pro- 
posed theory [191 ]. A related approach called interact- 
ing subsystem approximation (ISA) has been proposed in 
Ref. [32| for ASEPs with a single defect or several con- 
secutive defects (bottleneck) . Here it was suggested that 
due to the defect bonds there are 3 segments in the sys- 
tem: two homogeneous ASEPs are coupled by a segment 
that includes all sites that surround defect bonds. Ex- 
plicit results have been used inside the segments, and 
mean-field assumptions have been utilized for particle 
dynamics between the segments. A better agreement 
with Monte Carlo computer simulations has been found, 
and the method was also successfully applied to describe 
interactions of defects with boundaries. It was argued 
that ISA can be used for analyzing properties of general 
ASEPs with disorder [H,!!!]. However, ISA has not been 
applied for the systems with 2 defects at finite distances 
from each other, and because of this observation it is 
difficult to apply ISA for understanding mechanisms of 
more complex inhomogeneous asymmetric exclusion pro- 
cesses. A slightly different method of calculations has 
been proposed by Chou and Lakatos Q, who applied a 
finite segment mean-field theory (FSMFT). According to 
this approach, the segment of finite length n that cov- 
ers the defect and surrounding sites is considered, and 
its dynamics is fully described by solving explicitly for 
eigenvectors of the corresponding transition rate matrix. 
The segment is then coupled in the mean-field fashion to 
the rest of the system. However, this approach becomes 
numerically quite involved for cluster sizes larger than 
«20, and it also limits its applicability. Different stud- 
ies of asymmetric exclusion processes with disorder point 
out to importance of correlations in the system. It is rea- 
sonable to expect that correlations are stronger near the 
slow defect sites. However, it is not clear how far from 
the local inhomogeneity and how fast these correlations 
decay. In addition, it is also unclear how correlations 
from two close defects affect each other. The goal of this 
paper is to investigate the role of correlations in dynamics 
of ASEPs with disorder. By analyzing several analytical 
approaches in combination with extensive Monte Carlo 
computer simulations it will be shown that a successful 
description of disordered driven diffusive systems can be 
achieved by properly accounting for correlations near the 
defect bonds. 



II. MODEL AND THEORETICAL 
DESCRIPTION 

We investigate a totally asymmetric simple exclusion 
processes with disorder. In the one-dimensional lattice 
the particle at the site i can jump forward with the rate 
Pi if the next site i + 1 is unoccupied, otherwise it stays 



at the same place. The particle can enter the system 
with the rate a if the site is empty, and it can also exit 
the lattice with the rate /3. When all pi = 1 we have 
a homogeneous ASEP for which dynamic properties are 
known explicitly from exact solutions P, Q . ASEPs with 
disorder correspond to the situation when there is inho- 
mogeneities in the transition rates, and pi are drawn from 
arbitrary distributions. Numerous theoretical and com- 
putational studies indicate that in the limit of large times 
the dynamics in the system can be determined by com- 
paring entrance rate, exit rates and the transition rate 
at the slowest defect bond plj [33| . This observation has 
a significant consequence for properties of ASEPs with 
disorder, yielding a generic phase diagram with 3 phases. 
When the entrance is a rate-limiting process the system 
can be found in the low-density phase, while for slow ex- 
iting the high-density phase governs the system. If the 
rate-limiting process is the transition via the slowest de- 
fect bond the system is in the maximal-current phase. In 
this maximal-current phase, a segregation of density pro- 
file into macroscopic high and low regions occurs at the 
location of slowest defect bond. Other defects only per- 
turb the density profile on a local scale. However, when 
the number slowest defect bonds exceeds two or more the 
above picture needs modification. Furthermore, the pre- 
vious studies on disordered ASEPs lack investigations on 
correlation effects induced by defects. To address these 
questions, we analyze the simplest model with 2 identi- 
cal defects in the bulk of the system far away from the 
boundaries. It was shown earlier [32j that positioning 
of the slow defects close to the boundaries leads only to 
rescaling of the effective entrance and/or exit rates, and 
we will not consider this possibility in this paper. Note 
that in this paper we are using terms of defect bonds and 
defect sites. To clarify, we define the defect site as the 
site i from which the particle hopes to the site i + 1 with 
the rate q < 1. Correspondingly, the bond connecting 
sites i and i + 1 is a defect one. 



A. Defect Mean-Field Theory 

Consider a totally asymmetric exclusion processes with 
open boundaries and with 2 slow defective sites at i — d\ 
and i = d% at a distance d with d2 — d± = d (separated by 
d— 1 normal sites), as shown in Fig. 1. At the defects the 
particle jump to the right with the rate q < 1, in all other 
sites the hopping rate is equal to one. It can be seen that 
two defects divide the system into three segments. 

The particle dynamics inside each segment can be cal- 
culated exactly, however, it is assumed that there are no 
correlations between the segments. If entrance to the lat- 
tice is the slowest process then the system can be found 
in low-density (LD) phase with stationary current and 
bulk densities given by 

J = a(l-a), Pbuik = a. (1) 



3 



a 







• 


• • • 




• 












d 2 





FIG. 1: Fig.l: Schematic of ASEP with two defective sites 
at i = di and i — d2 separated by d — 1 normal sites i.e., 
di — d\ = d. The reduced hopping rates at each defect is 
equal to q. In normal sites the hoping rates are one. 



Similarly, when the exit becomes a bottleneck process 
the system is in high-density (HD) phase with 



J = /3(l-p), Pbu i k = l-p. 



(2) 



Note that in both phases particle densities near the de- 
fect bonds will deviate from the bulk values. The more 
interesting case is when the dynamics in the system is 
governed by transitions via local inhomogcneitics. In 
this phase, which has the maximal current, we expect 
to have density phase segregation similar to the case a 
single defect ASEPs. We emphasize that all our inves- 
tigation in this paper is on this maximal-current phase. 
First consider a lattice segment after the second defect. 
The dynamics in this part of the system is controlled 
by the entrance of particle via the defect, then it has a 
low-density profile with unknown bulk density p* < 1/2. 
Similar arguments can be used to analyze the density pro- 
file in the segment before the first defect. Here the flux 
is limited by the exit rate via the local inhomogeneity, 
leading to the high-density phase. Since at stationary- 
state condition the flux through any segment should be 
the same, J = p*(l — p*), the bulk density in this seg- 
ment is equal to 1 — p* . The region between two defects 
can be viewed as asymmetric exclusion process on a fi- 
nite lattice with d sites. The effective entrance and exit 
rates to this segment can be easily evaluated using our 
mean- field assumptions, 



oteff = Peff =q(l-p*). 



(3) 



The stationary properties of the lattice segment with 
d sites between the defects can be evaluated explicitly by 
utilizing exact results for finite-size ASEPs [5|. Specifi- 
cally, the particle flux is given by 



T( a* Rd-i(l/P)-Rd-i(l/a) 
J ° (a <^ R d ( m -R d (l/a) ' (4) 

where the function Rd{x) is defined as 
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H d{x)-2_^ d u d+1 _ ), x ■ 



(5) 



To understand the density profile in the segment we 
can use a domain-wall picture of asymmetric exclusion 
processes [lij . Since the entrance and exit rates are the 
same [see Eq. ©], the domain wall that separates high- 
density and low-density blocks performs an unbiased ran- 
dom, leading to a linear density profile with a positive 
slope. Explicit expressions for particle densities can also 
be found in Ref. [5[. The full description of dynamics 
in ASEPs with two defects is obtained by solving for the 
unknown parameter p*. It can be done by applying the 
condition of stationarity in the particle flux, 



J = p*(l - p*) = J (a e f f ,P ef f,d). 



(6) 



This equation can always be solved analytically or nu- 
merically exactly for any number of sites between local 
inhomogeneities, leading to stationary particle currents 
and density profiles, it is important to note that there is 
a particle-hole symmetry in the system because defects 
are far away from the boundaries. To illustrate our ap- 
proach let us calculate dynamic properties of ASEPs with 
two defects for several values of the parameter d. First, 
let us analyze the simplest case of d = 1 with consecutive 
defects in the bulk. It can be shown that for this system 



J (a,/3,d = 1) 



a/3 
a + /3' 



Then Eq. ([6|) can be written as 



p*(l-p*)=q(l-p*)/2, 



(7) 



(8) 



which produces simple expressions for the bulk density 
and the particle current, 



p*=q/2, J = q(2-q)/i. 



(9) 



The density I at the site between the defects can also 
be found from the condition that the flux via this site, 
J = ql(l — p*), should be equal to the flux through other 
segments, and this yields 



I = p*/q= 1/2. 



(10) 



This result could also be obtained from the particle- 
hole symmetry arguments. Note that for q = 1 we ob- 
tain p* = 1/2 and J = 1/4 as expected for homogeneous 
ASEPs in the maximal-current phase. For d = 2 there 
are two lattice sites between the defects, and stationary 
properties of this system can also be obtained analyti- 
cally. From Eq. ([5]) one can easily derive 



R 2 (x) 



(11) 



which produces the following expression for the current 
in the segment between the defects 
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J {a,(3,d = 2) 



i + i 



1 
a/3 



(12) 



Using the expression for the effective entrance and exit 
rates for the segment between the inhomogeneities [see 
Eq. ©], the condition for the stationary current leads 
to 



P*(l-P*) = 7 



2q(l-p* 



3 + 2q(l-p*) 



(13) 



This quadratic equation can be solved, and taking the 
physically reasonable root we obtain 



2q + 3 - ^9 + 12q- I2q 2 



(14) 



T _ 8g 2 - 6g - 9 + 3 ^9 + \2g - 12g 2 

It can be checked that for q = 1 these equations reduce 
to expected relations p* = 1/2 and J = 1/4. We can also 
calculate the densities l\ and I2 at the sites between the 
defects. Because of the particle-hole symmetry one can 
argue that 



(16) 



and the density at the first site can be found by ana- 
lyzing the current via the first defect, 



J = q(l-p*)(l-h) = p*(l-p*). (17) 
Then we have 



h = 1 



f - 2q - 3 + V9 + I2q- I2q 2 



(18) 



We have solved equation (6) for the case d = 3. In this 
case we have: 



R 3 (x) = 2x 2 + 2x 3 + x 4 



(19) 



After some lengthy but straightforward algebra we ar- 
rive at the following cubic equation for p*\ 



4q A (p*y-(8q' + 6q){p*y + (6q 2 + 6q + 4:)p*-q(3 + 2q) = 0. 

For brevity we avoid writing the answer explicitly. An- 
alytical results for ASEP with 2 defects can also be ob- 
tained in the limit of very large distances between the 




FIG. 2: Fig.2: (Colour online) J vs q for d = 1, 2, 3 obtained 
by DMF method and MC simulation. 



inhomogeneities (d 3> 1). In this case the segment be- 
tween the defects can be viewed as a homogeneous ASEP 
in the state of the phase transition between high-density 
and low-density phases (a e ff = Peff)- This corresponds 
to a linear density profile for the segment between the 
defects. Then it leads to the following expression for the 
current 



p*0-p*) = q(l-p*)[l-q(l-p*)], 



and finally we obtain 



J = 



(l + q) 



2 ' 



(21) 



(22) 



These results are identical to stationary properties of 
ASEP with only one local inhomogeneity far away from 
the boundaries obtained using DMF approximation [l|| , 
suggesting that 2 defects at large distances do not affect 
each other [8j. For a general d equation (6) leads to a 
polynomial equation of order d for the unknown p* . For 
d > 3 this equation can be solved numerically to find the 
acceptable answer. In figure (2) we have sketched the 
behavior of current J as a function of q for d = 1, 2 and 
3 and have compared them to the results obtained via 
Monte Carlo simulations. As expected J is an increasing 
function of both q and d. DMF notably underestimates 
the current in comparison to the MC simulation espe- 
cially in the intermediate values of q. 



B. Interacting Subsystem Approximation 

Interacting subsystem approximation (ISA) is another 
method of calculating stationary properties of ASEPs 
with a single defect or a single bottleneck developed by 
Greulich and Schadschneider [33|. It can be easily ex- 
tended to the case of asymmetric exclusion processes with 
2 defects separated by d lattice sites. Similarly to DMF 
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J12 = (l-/i)(l-^*)=a 2 (l-p*), 



(25) 



where a 2 is the effective rate to enter the second seg- 
ment. At large times we expect to find the system in the 
stationary state, i.e., J — J12, yielding 



FIG. 3: Fig. 2: (Colour online) Fig. 3: (Colour online) inter- 
acting subsystems connected via mean-field assumption. For 
d>2 the system is divided into five segments. 



this method divides the lattice into several segments. 
Particle dynamics inside the segments is treated exactly, 
while between the segments mean-field assumptions are 
made. ISA differs from DMF in the defining of segments. 
In DMF the position of defects separates different parts, 
and there is always 3 segments in the system. In ISA 
the sites that are connected by the defect bond are put 
together in one segment, as shown in Fig. 3. 

For d = 1 there are also 3 parts in the lattice, and 
the middle segment has 3 sites. For d = 2 there are 
4 segments and 2 middle segments (with 2 lattice sites 
each) border each other. For any larger distance between 
local inhomogeneities ISA assumes 5 segments: see Fig. 
3. Note that the size of the middle segment is equal to 
d-2. 

Let us consider a general case of 5 segments (d > 2) 
for computation of stationary properties of ASEPs with 
2 defects. As was argued above, the system can be found 
in one of three phases: LD, HD or HD/LD (maximal- 
current). Since the derivation of properties in HD and 
LD phases is the same as for DMF approach, we concen- 
trate on description of the maximal-current phase. As 
before we assume that the bulk density in the segments 
1 and 5 are 1 — p* and p* correspondingly. Let us define 
li and I2 as the probabilities to find the particles at the 
corresponding sites of the segment around the first de- 
fect. Similarly, I3 and I4 describe densities in the segment 
around the second defect bond. For the middle segment 
with d — 2 lattice sites we define Xi for % = 1, • • • , d — 2 
as the particle density at i-th site of this segment. As for 
DMF approach, the existing particle-hole symmetry sim- 
plifies calculations significantly. Specifically, it suggests 
that 



h = l-lx, h = 1-/2, Xi = 1 - Xd+x-i- (23) 

The overall particle current in the system can be writ- 
ten as 



J = p*(l-p*), (24) 

while due to the mean-field assumptions the current 
between the first and the second segments is equal to 



1 - a 2 = h = (1 - P*)- (26) 

The current between segments 2 and 3 can be pre- 
sented in the several ways, 

J 23 - Z 2 (l ~ x x ) = foh = a 3 (l - Xl ), (27) 

with 82 being the effective exit rate from the segment 
2, while «3 is the effective rate to enter the segment 3. 
When the system reaches stationary phase, J = J23, and 
we obtain 

P2 = 1-2:1, a 3 = — . (28) 

1 — Xl 

Because of the particle-hole symmetry the effective en- 
trance and exit rates from the segment 3 are the same, 
<X3 = /?3- The particle current via the ASEP segment 
with TV sites and with entrance and exit rates a and 8, 
respectively, J(a, 3, N), can be calculated explicitly [a]. 
Then to obtain stationary properties of ASEP with 2 de- 
fects in the maximal-current phase the following system 
of equations should be solved, 

,1 — P* 1 — Xi 
p*(l-p*) = qj( --, i,2); 



m-p*) = J{^^,^^,d-2). (29) 
1 — X\ 1 — Xl 

where p* and Xx are 2 unknown variables. The expres- 
sion on the right side of the first equation describes the 
current inside the segment 2 and 4. Because the hopping 
rate is q < 1, the effective entrance and exit rates must 
be rescaled by the same factor. The right side of the sec- 
ond equation gives the current inside the segment 3. The 
application of ISA for d — 1 and d = 2 cases is different. 
In the case of 2 consecutive defect bonds the system is 
divided only in 3 segments. The middle segment has 3 
sites that surround defect bonds. In the HD/LD phase 
the effective entrance rate is a 2 = 1 — p*, and the sta- 
tionary properties can be obtained by solving only one 
equation 

p*(l-p*)= g j(l^-i,i-^l,3). (30) 

q q 
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Using Eqs. (|4]) and ([5]) for the middle segment with 
equal entrance and exit rates gives us 



p*(l-p*) = 



q(l - p*) [2(1 - p*) + 3q] 



(31) 



2[2(l-p*)2 + 3( 7 (l-p*) + 2( Z 2] 
which can be simplified into the following expression, 
A(p*) 3 -2(3q + A)(p*) 2 +A(q + l) 2 p*-q(3q + 2) = 0. (32) 
This cubic equation can be solved explicitly, yielding 
p* = [3q + 4 + (4 - 3q 2 )/D + D] /6, (33) 

where 



D = 



-8 + 9q 2 - 27q 3 + 3-y/3 Vl6t? 3 - <? 4 - 18<? 5 + 28<? £ 



(34) 

ISA also works differently in the case of d = 2. There 
are 4 segments in the system, and because of the neglect 
of correlations between segments 2 and 3 we have 




ISA, d = 1 
ISA, d = 2 
ISA, d = large 
MC, d = 1 
MC, d = 2 
MF, d = large 



0.25 0.5 0.75 1 

1/3 q 

FIG. 4: Fig. 4: (Color online) Current vs q for d — 1, 2 and a 
large d obtained within ISA method and MC simulation. In 
simulations we have taken a = ft = 0.6. 



l 2 (l-l 3 ) = l 2 2 = p*(l-p*). 



(35) 



Then the effective entrance rate to the segment 2 is 
a% = 1 — p*. and the effective exit rate is equal to 
(3-2 = h — y/p*0-~ P*)- The unknown parameter p* is 
determined from the equation for the stationary current, 



(36) 



Substituting the values of the effective boundary rates 
and utilizing Eq. (|12[) we obtain 



III. CORRELATIONS NEAR DEFECT 
A. Monte Carlo Simulations 

In this section we aim to investigate correlations in 
the vicinity of defects. We restrict ourselves to adjacent 
two-point correlations and will present our results for the 
general two-point and multi-point correlations in a future 
work. Let us now introduce the normalized connected 
two-point correlation function d between the neighbour- 
ing sites i and i + 1. This quantity is defined as follows: 



p*^T^7+(p*f /2 =qV^7- (37) 
which can be recast in the form of a cubic equation: 



2(p*r - (2q + l)(p*Y + (q 2 + 2q)p* - q 2 = 0, (38) 

This equation can be solved analytically but for brevity 
we do not write the solution. It can also be shown that in 
the limit of d ^> 1 ISA method with 2 defects produces 
the stationary current and bulk densities which are in- 
distinguishable form the situation with only one defect 

H, 



3q + 2 - V V - 4q + 4 /4 



J 



2q - 9q 2 



Aq 



/8. (39) 



Let us now exhibit the dependence of J on q in the ISA 
method. In figure (4) we have drawn J vs q for d = 1, 2 
and have compared the results to those obtained by DMF 
method and MC simulations. 

In general, ISA method gives a better estimation of 
current compared to DMF at least for small values of d 
we have considered. 



Ci 



H+l) 



i=l,... ,L-X&0) 



i+i/ 



H+l) 



The function Ci measures the correlation and it lies 
between —1 and 1. Negative values correspond to anti- 
correlation between neighboring sites whereas a positive 
value signifies correlation. The values near zero are re- 
garded as uncorrelated. Fig. (5) depicts the simulated 
profiles of correlation at d — 10 and 100 each for three 
values of q. The system size is L = 500 and we have 
taken a = (3 = 0.6 in all our simulation results unless 
stated otherwise. The system has been updated for T 
Monte Carlo steps. Each step consists of L moves. In 
each move, we randomly choose a site and update its 
status according to ASEP rules described above. We dis- 
card the first steps to ensure reaching steady state, and 
we have accumulated data separated by 10 MC steps to 
avoid any possible temporal correlations. The value of T 
is taken 10 s in our simulations. 

Two defects are symmetrically placed with respect to 
chain mid point. We observed that correlations are large 
in sites between the defects. There is a rather strong 
anti-correlation in the sites immediately after the first 
defect and before the second defect. The correlations are 
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site number site number 



FIG. 5: Fig. 5: (Color online) Profile of normalized correlation 
function at d = 10 (top) and d = 100 (bottom) for q — 
0.1,0.3,0.5. 



FIG. 6: Fig.6: (Color online) Profile of normalized correlation 
functions at various values of d for q = 0.1 (top) and q = 0.3 
(bottom). 



growing up for middle sites where the maximum value 
is achieved. It can be seen that correlations are greater 
when d is increased. This is unexpected and counterintu- 
itive because increasing the distance between the defects 
reduces their interaction. It has been observed via MC 
simulations that when d is increased the current reaches 
asymptotically to its mean- field value Jmf — rf+^p 
[I3,[3l[. Therefore, one expects the correlations to ex- 
hibit a reducing behavior with respect to distance d but 
this is not observed in our simulations. To have a deeper 
understanding, we have sketched the behavior of correla- 
tion profile upon varying the distance d for q = 0.1 and 
q = 0.3 in Fig. 6. 

For fixed values of q, increasing the distance d be- 
tween the defects gives rise to enhancement of corre- 
lations/antocorrelations. For instance, the correlation 
value in the middle point rises from roughly 0.5 at small 
d ~ 10 to 0.65 for d ~ 100. It can be observed that 
there is no notable difference in correlation values for 
d larger than 100. Moreover, the correlations are al- 
ways greater than anti-correlations. By increasing q, the 
correlations/anti-correlations are notably reduced in val- 
ues. This is expected since in the limit of homogeneous 
ASEP where q — > 1 the correlation functions become very 



small. Here we wish to make a pause and have a discus- 
sion on correlations in normal ASEPs. In fact the middle 
segment between two defects can be regarded as an ASEP 
chain with length d with equal entrance and exit rates. 
To the best of our knowledge, correlations in ASEP with 
random sequential update, has only been analytically dis- 
cussed by Derrida and Evans who obtained exact analyt- 
ical expression for a general two-point function and made 
a conjecture to generalize their findings to n-point func- 
tion [39(. Their study was restricted to the special case 
a — (3 = 1 and they found that long range correlations 
persist in the bulk which was attributed as a boundary 
effect. In order to see if the large value of the connected 
two-point function survives in the normal ASEP with 
equal entrance and exit rates, we performed MC simula- 
tions. Our results show that when a = (3, the profile of 
Ci reaches a small constant (almost zero) in the bulk. 
The correlations become large near boundaries. This 
boundary behaviour depends on whether a = (3 < 0.5 
or a — (3 > 0.5. Figure (7) illustrates this aspect. 

We recall that correlations in other types of update 
such as parallel updating has been discussed in [3, [TtJ ■ 
It is worthwhile to examine the behavior of density pro- 
file between defects. Dong et al have recently shown via 
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FIG. 7: Fig. 7: (Color online) Profile of normalized correlation 
functions in a normal ASEP chain with a = (3 . 
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FIG. 9: Fig. 9: (Color online) Dependence of normalized cor- 
relation functions at three selected sites on q (top) and on d 
(bottom). 



FIG. 8: Fig. 8: (Color online) Profile of density at various 
values of d for q — 0.1. The profile exhibits a linear behaviour 
with positive slope. This behaviour is associated to the exis- 
tence of wandering shock in the region between left and right 
defects. 



extensive MC simulations that the density profile takes a 
linear shape between defects [T3| ■ This behavior remains 
unchanged in ASEP with extended objects [3l|. For the 
sake of completeness, we show some typical density pro- 
files in Fig. 8. 

The interesting point is the absence of boundary layer 
in this phase-segregated regime. It would be illustrative 
to look at the dependence of two-point correlation func- 
tions at some particular sites on values of q and d. These 
results are sketched in Fig. 9 where correlations at the 
first defect site (d\), its rightmost site {d\ +1) and in the 
middle site of the chain are considered. 

Note that all correlations/anti-correlations approach 
zero when q tends to one. Moreover, increasing defect's 
separation d increases the correlations. The dependence 
on d is more interesting. While the values of correlation 
functions reach an asymptotic value at large d, the behav- 
ior is not monotonous. Especially for C^+i correlation 
increases up to a maximum and then begin to decrease 



smoothly towards its asymptotic value. The value of d 
where C^+i is maximized does not show a significant 
dependence on q. In can be concluded that varying d 
can dramatically affect the system characteristics as far 
as correlations are considered. 



B. Analytical theory 

Our simulation findings in the preceding section con- 
firms that in between the defects the correlations are no- 
tably higher than other sites. In this section we try to 
develop a theoretical framework to capture this feature. 
Suppose we have two slow defects located in the bulk at 
sites k and I (k < I) respectively both with rate q. The 
mean occupation at site i in the steady state is denoted 
by n i = in which n = 0, 1 is the occupation number 
at site i. We assume that a simple MF assumption, i.e., 
(7"tTt+i) = (Tj)(Tj-|_i) = riini + i holds for all sites except 
i = k, ■ ■ ■ ,1 i.e.; defective sites themselves and all the 
sites between them. At these sites the correlations are 
strong enough to violate the simple mean-field assump- 
tion. Let us introduce two-point correlation functions rrii 
in the following way, 
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HH+lj 



(41) 



There are L 



I - k 



2 unknowns, namely, 
ni,n2,--- , til, rrifc, TOfc+i, • • • , m; and lastly the current 
J. In the stationary state there exists L + 1 equations 
among these unknowns. Let us label them by Ao to Al ■ 
These equations can be obtained by expressing the cur- 
rent J in terms of the mean site densities and two point 
correlators. The first equation, Aq, is J = — ni). The 
equations Ai for i = 1, • • • , k — 1 and i = I + 1 , • • • ,L—1 
have the following form, 



J = r H {l-n i+1 ). (42) 

Equation Afe is: 

J = q(r k (l - r k+1 )) = q(n k - m k ). (43) 
Similarly, equation A; is given below 

J = q(n(l - n +1 )) =q{m-mi). (44) 



Equations Ai (i = k + 1, 
form 



, / — 1 ) have the following 



J = (Ti(l - r i+ i)) = (n, - m i+ i). (45) 
lastly equation Al is as follows, 

J = /3(r L > = (3n L . (46) 

We do not intend to add more equations. Unfortu- 
nately the above equations are nonlinear and it would 
be a formidable task to solve them analytically. Alterna- 
tively, we shall utilize a numerical approach to solve the 
system of equations by exploiting their recursive struc- 
ture. This approach was originally introduced in 29] in 
the context of disordered ASEP and was later applied to 
the problem of two intersecting ASEP chains [30]. Ac- 
cording to this numerical scheme, we assign a value to J. 
Having </, it is possible to iterate equations (in forward 
direction) and obtain n\ up to n k . Then we proceed to 
find m k by incorporating the relation J = q(n k — m k ). 
At this stage it is not possible to proceed further because 
both n k+ i and m k+ i are unknown and we have only one 
relation between them : J = n k+ \ — m k +\. In order to 
proceed, we approximate n k+ i in the following way. Con- 
sider the rate equation for the 2-point function (r k T k+ i) 
which is governed by the following master equation: 



d(T k T k+1 ) 

dt 



(r fe _i(l - T k )r k+ i) - (T k r k+ i(l - T fe+2 )^.47) 



In the steady state, the left hand side becomes zero, 
and therefore two terms on the right hand side will be 



equal. To proceed further we have to approximate 3- 
point functions. This is achieved by utilizing the cluster 
mean- field assumption [40]. According to this assump- 
tion we replace any three point function by the product 
of 2-point functions as follows, 



{nin 3 n k ) 



(ninj}{njrik) 



(48) 



We then replace all 2-point functions by the product 
of 1-point functions except m k — (T k T k+ i). Then it is 
possible to express 1 — n k+ 2 in terms of n k ^i,m,i and 
n k+ i. After some algebra a quadratic equation for n k+ i 
is obtained: 



n k -\n 2 k+1 - m k n k ^in k+ i - m k J = 0. 



(49) 



The physically reasonable solution for this equation is 
given by 



m k 



n k+ i 



mi + 



im k J 



(50) 



Now it is possible to find mfc+i via equation J = n k+ i — 
m k +\. Analogous to the above procedure we can find 
n k+ 2 as follows: 



n k +2 = 



qm k m k+1 + yj (qm k m k+ i) 2 + AqJn k n 2 k+1 m k+ _ 



2qn k n k+1 



{51) 



After having n k+ 2 we simply obtain m k+ 2 via equation 
J = n k+ 2 — fn k+ 2- Now it would be possible to proceed 
iteratively after taking into account some approximation. 
To this end we recall the equality 



d(nTi- 



dt 



= (n-xii - n) n+1 ) - (n n+1 (i - r i+2 )).(52) 



Putting the left hand side equal to zero in the steady 
state, utilizing cluster mean-field in three point functions 
and finally substituting rrii + i by rrii+i — n,; + i — J we 
arrive at the following equation: 



H+l — 



2n l - 1 n i 



(53) 



Note that we have approximated (Ti_iT J+ i) by the 
mean-field relation (t^-i) (ri + i). We can iterate equation 
(53) together with nii+i = rii + i — J from i = k + 2 to 
I — 2 to find the corresponding ni and up to i — I — 1 . 
The site i = I needs to be treated separately. Following 
the same strategy we easily find: 



m 



{54) 
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FIG. 10: Fig. 10: (Color online) Current vs q for d = 2 
obtained by various analytical and numerical methods and 
MC simulation. J approaches to 0.25 when q tends to one in 
accordance to normal ASEP. 



FIG. 11: Fig.ll: (Color online) Current vs d for q = 0.1, 0.3 
and 0.5 obtained by the numeric method and MC simulation. 



From which one can compute mi. In a similar fashion, 
we can obtain nz+i. We only should shift up all the sub- 
scripts in equation (54) by one. Now it is possible again 
to proceed iteratively to the end of the chain and evaluate 
riL which gives us the output current J out . If the given 
value of input J were correct, the output current J out , 
which is [3til, should be the same, up to a given precision, 
as the input value of J. By systematically increasing the 
input J in an small amount 5 J, we can determine the cor- 
rect J and correspondingly the mean densities n\ , • • • ,til 
together with correlators mk, ■ ■ ■ ,mi. In Fig. (10) the 
dependence J on q obtained by the numeric scheme de- 
vised above is sketched. For the sake of comparison, we 
have augmented the figure with the analogous graphs ob- 
tained by MF, DMF, ISA and MC methods. 

The result of the numerical scheme is almost identical 
to ISA method. They both are in very good agreement 
with Monte Carlo simulations. However, the numerical 
scheme has an advantage over ISA method in the sense 
that it can easily be implemented for any d, whereas find- 
ing the solution of the ISA nonlinear set of equations, 
i.e., Eq. (29) is not an easy task even by employing 
advanced numerical methods. Finally in figure (11) we 
have sketched the dependence of J on d obtained from 
MC simulation and the numeric scheme. 

The results of the numeric method are in rather good 
agreement with those obtained by MC simulations. J is 
an increasing function of d and becomes saturated after 
some short g-dependent distance. The results confirms 
the earlier finding in [lOj. Note that the length scale on 
which J recovers its single-defect value is of the same or- 
der of magnitude of the correlation length in the density 
profile near boundaries. Finally we would like to add 
that our numerical scheme is not capable of reproducing 
the profile of correlators obtained via MC simulations. 
Figure (12) depicts the profile of unnormalised adjacent 
two-point correlation function (rj+iTj) — (Tj_|_i)(Tj) for two 
methods of simulation and numerical scheme. 
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FIG. 12: Fig. 12: (Color online) Profile of corrrelators for 
d = 100 and q = 0.3. 



We see that within the numerical framework both the 
value and the extension of correlators are small in com- 
parison to the simulation results. The reason lies in the 
implementation of a series of approximation in several 
places in this numerical algorithm. 



IV. SUMMARY AND CONCLUSIONS 

An open ASEP chain with two defective sites with re- 
duced hopping rates q < 1 has been investigated. The 
system current and mean site densities at defective sites 
and their vicinities have been obtained by two analyti- 
cal methods namely defect mean-field (DMF) and inter- 
acting subsystem approximation (ISA). Both methods 
combine mean-field approach near defects with known 
exact solutions. Our results are accomplished by exten- 
sive Monte Carlo simulations. We focus on the phase- 
segregated phase in which defects globally affect the sys- 
tem properties and the system is not input/output rate- 
limited but rather defect-limited. MC simulations have 
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revealed strong short ranged correlations at the defec- 
tive sites and at all the site between them. Additionally, 
the profile of density between defects takes a linear form 
which marks the existence wandering shock in this inter- 
mediate region. In order to take into account these corre- 
lations, we have introduced a numerical approach which 
utilizes a cluster mean-field assumption. Comparison of 
the three methods show that ISA and the numeric ap- 
proach give a current value which is in a good agreement 
with MC simulations. DMF method however, only gives 
a good results compared to MC and other two methods 
for small q less than 0.1 which is due to strong corre- 
lations. Furthermore, we have obtained the profile of 
neighbouring two-point correlation function throughout 
the chain. It is shown that these two point correlators 
exhibit a rather strong anti correlation at the first defec- 
tive site then they grow to the middle of the defects and 
after that they start diminishing. In general, two point 
correlators will tend to a tiny value when q approaches 
one. However, the dependence of two point correlators 



for a fixed q as a function of distance between defects 
is not monotonous. Despite reaching to an asymptotic 
value for large distances, one observes a peak at short 
distances for the two point correlators near defects. Our 
theoretical analysis indicates that correlations are criti- 
cally important for dynamics of particles in disordered 
ASEPs. It also shows that it is possible to devise an ap- 
proximate method that can take into account these corre- 
lations, providing a satisfactory description of stationary 
properties. 
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